Remove contaminants and filter samples with less than 1000 reads

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../source/main_packages.R")
source("../../source/functions.R")

# Load supplementary packages
packages <- c("decontam", "kableExtra", "microbiome", "cowplot")
invisible(lapply(packages, require, character.only = TRUE))

# create output folders if needed
if (!dir.exists("../../../output/1_MED/1D")) {dir.create("../../../output/1_MED/1D")}

Import MED phyloseq object

ps <- readRDS("../../../output/1_MED/1C/1C_MED_phyloseq.rds")

Run Decontam

Detect contaminants

# set threshold
threshold=0.6
  
# Preprocess
df <- as.data.frame(sample_data(ps)) # Put sample_data into a ggplot-friendly data.frame
df <- data.frame(sample_data(ps)) # Put sample_data into a ggplot-friendly data.frame
df <- df[order(df$Read_depth),]
df$Index <- seq(nrow(df))

# read depth 
p <- ggplot(data=df, aes(x=Index, y=Read_depth, color=Control)) + geom_point()
p

# Identify contaminants - Prevalence method
sample_data(ps)$is.neg <- sample_data(ps)$Control == "Control sample"
contamdf.prev <- isContaminant(ps, method="prevalence", neg="is.neg", threshold=threshold)

# table of positive and false contaminant
table(contamdf.prev$contaminant)

## 
## FALSE  TRUE 
##    67     3

# contaminant nodes
decontam_asv_MED  <- row.names(contamdf.prev[contamdf.prev$contaminant == TRUE, ])

# make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(ps, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$Control == "Control sample", ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$Control == "True sample", ps.pa)

# make dataframe from phyloseq objects of presence-absence
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev$contaminant)

# plot
p2 <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
  xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")

p2

# write the contaminant ASV table detected with prevalence
tax <- tax_table(ps) %>% as.matrix()
decontam_df_MED <- tax[row.names(tax) %in% decontam_asv_MED , ]

# print contaminants
decontam_df_MED   %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Kingdom	Phylum	Class	Order	Family	Genus	Species
N0005	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Moraxellaceae	Enhydrobacter	NA
N0315	Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacterales	Yersiniaceae	Rahnella1	NA
N0852	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA

# histogram of the scores
p <- hist(contamdf.prev$p, 100, ylim = c(0,25), xlim = c(0,1), main="", xlab="p", ylab="Frequency")

Save contaminant table

Remove contaminants from the phyloseq object

# remove contaminants ASV
alltaxa <- taxa_names(ps)
decontam_taxa <- alltaxa[!(alltaxa %in% decontam_asv_MED)]
ps2 <- prune_taxa(decontam_taxa, ps)

# check ps 
ps2 <- check_ps(ps2)
ps2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 67 taxa and 123 samples ]
## sample_data() Sample Data:       [ 123 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 67 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 67 reference sequences ]

compute_read_counts(ps2)

## [1] 6469784

compute_read_counts(ps)-compute_read_counts(ps2)

## [1] 17460

ps <- ps2

Remove blanks

# blanks
ps.blanks <- subset_samples(ps, Strain=="Blank")
ps.blanks <- check_ps(ps.blanks)
ps.blanks

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 24 taxa and 8 samples ]
## sample_data() Sample Data:       [ 8 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 24 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 24 reference sequences ]

compute_read_counts(ps.blanks)

## [1] 16211

# supprimer blanks
ps.filter <- subset_samples(ps, Strain!="Blank")
ps.filter <- check_ps(ps.filter)
ps.filter

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 67 taxa and 115 samples ]
## sample_data() Sample Data:       [ 115 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 67 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 67 reference sequences ]

compute_read_counts(ps.filter)

## [1] 6453573

# check nreads
summarize_phyloseq(ps.filter)

## Compositional = NO2

## 1] Min. number of reads = 952] Max. number of reads = 7361593] Total number of reads = 64535734] Average number of reads = 56118.02608695655] Median number of reads = 251717] Sparsity = 0.6486696950032456] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons 
##         (i.e. exactly one read detected across all samples)010] Number of sample variables are: 12SampleStrainFieldCountryOrganSpeciesRunControlSpecies_italicStrain_italicRead_depthis.neg2

## [[1]]
## [1] "1] Min. number of reads = 95"
## 
## [[2]]
## [1] "2] Max. number of reads = 736159"
## 
## [[3]]
## [1] "3] Total number of reads = 6453573"
## 
## [[4]]
## [1] "4] Average number of reads = 56118.0260869565"
## 
## [[5]]
## [1] "5] Median number of reads = 25171"
## 
## [[6]]
## [1] "7] Sparsity = 0.648669695003245"
## 
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
## 
## [[8]]
## [1] "8] Number of singletons = 0"
## 
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n        (i.e. exactly one read detected across all samples)0"
## 
## [[10]]
## [1] "10] Number of sample variables are: 12"
## 
## [[11]]
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"

Counts

Create dataframe

x <- c("Culex pipiens", "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar", 
       "Culex quinquefasciatus", "Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)",
       "Aedes aegypti (Guadeloupe)",
       "Total")
y <- c("Reads", "Oligotypes", "Samples")

df <- data.frame(matrix(ncol=3, nrow=9))
rownames(df) <- x
colnames(df) <- y

df2 <- df
df3 <- df

Whole + ovary

# Culex pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)

## All Strain
### oligotype
nrow(ps.pipiens@otu_table) -> df["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data) -> df["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens) -> df["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Reads"]


# Culex quinquefasciatus
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)

## All Strain
### oligotype
nrow(ps.quinque@otu_table) -> df["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data) -> df["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque) -> df["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque %>% 
                 subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
                 check_ps())) -> df["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Reads"]


# Aedes aegypti
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Reads"]

# Total
### oligotype
nrow(ps.filter@otu_table) -> df["Total", "Oligotypes"]
### samples
nrow(ps.filter@sam_data) -> df["Total", "Samples"]
### reads
compute_read_counts(ps.filter) -> df["Total", "Reads"]

df %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	2574603	49	83
Field - Bosc	803574	36	23
Field - Camping Europe	869617	39	25
Laboratory - Lavar	901412	48	35
Culex quinquefasciatus	1925054	59	21
Field - Guadeloupe	1721053	47	7
Laboratory - Slab TC (Wolbachia -)	204001	45	14
Aedes aegypti (Guadeloupe)	1953916	54	11
Total	6453573	67	115

Whole

# Culex pipiens
ps.pipiens.whole <- ps.pipiens %>% subset_samples(Organ=="Whole")
ps.pipiens.whole <- ps.pipiens.whole %>% check_ps()

## All Strain
### oligotype
nrow(ps.pipiens.whole@otu_table) -> df2["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data) -> df2["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole) -> df2["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df2["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df2["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df2["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Reads"]



# Culex quinquefasciatus
ps.quinque.whole <- subset_samples(ps.quinque, Organ=="Whole")
ps.quinque.whole <- check_ps(ps.quinque.whole)

## All Strain
### oligotype
nrow(ps.quinque.whole@otu_table) -> df2["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data) -> df2["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.whole) -> df2["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.whole %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df2["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque.whole %>% 
                 subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
                 check_ps())) -> df2["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Reads"]



# Aedes aegypti
ps.aedes.whole <- subset_samples(ps.aedes, Organ=="Whole")
ps.aedes.whole <- check_ps(ps.aedes.whole)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.whole %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df2["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Reads"]


# Total
ps.whole <- ps.filter %>% subset_samples(Organ=="Whole")
ps.whole <- ps.whole %>% check_ps()
### oligotype
nrow(ps.whole@otu_table) -> df2["Total", "Oligotypes"]
### samples
nrow(ps.whole@sam_data) -> df2["Total", "Samples"]
### reads
compute_read_counts(ps.whole) -> df2["Total", "Reads"]

df2 %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	1251766	45	42
Field - Bosc	545790	33	13
Field - Camping Europe	98804	34	7
Laboratory - Lavar	607172	45	22
Culex quinquefasciatus	205761	57	18
Field - Guadeloupe	1760	38	4
Laboratory - Slab TC (Wolbachia -)	204001	45	14
Aedes aegypti (Guadeloupe)	1945807	54	9
Total	3403334	67	69

Ovary

# Culex pipiens
ps.pipiens.ovary <- ps.pipiens %>% subset_samples(Organ=="Ovary")
ps.pipiens.ovary <- ps.pipiens.ovary %>% check_ps()

## All Strain
### oligotype
nrow(ps.pipiens.ovary@otu_table) -> df3["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data) -> df3["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary) -> df3["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df3["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df3["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df3["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Reads"]



# Culex quinquefasciatus
ps.quinque.ovary <- subset_samples(ps.quinque, Organ=="Ovary")
ps.quinque.ovary <- check_ps(ps.quinque.ovary)

## All Strain
### oligotype
nrow(ps.quinque.ovary@otu_table) -> df3["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data) -> df3["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary) -> df3["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.ovary %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df3["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
# nrow(otu_table(ps.quinque.ovary %>% 
#                  subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
#                  check_ps())) -> df3["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
# nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Samples"]
# ### reads
# compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Reads"]



# Aedes aegypti
ps.aedes.ovary <- subset_samples(ps.aedes, Organ=="Ovary")
ps.aedes.ovary <- check_ps(ps.aedes.ovary)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.ovary %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df3["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Reads"]


# Total
ps.ovary <- ps.filter %>% subset_samples(Organ=="Ovary")
ps.ovary <- ps.ovary %>% check_ps()
### oligotype
nrow(ps.ovary@otu_table) -> df3["Total", "Oligotypes"]
### samples
nrow(ps.ovary@sam_data) -> df3["Total", "Samples"]
### reads
compute_read_counts(ps.ovary) -> df3["Total", "Reads"]

df3[is.na(df3)] <- 0

df3 %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	1322837	34	41
Field - Bosc	257784	26	10
Field - Camping Europe	770813	27	18
Laboratory - Lavar	294240	32	13
Culex quinquefasciatus	1719293	26	3
Field - Guadeloupe	1719293	26	3
Laboratory - Slab TC (Wolbachia -)	0	0	0
Aedes aegypti (Guadeloupe)	8109	30	2
Total	3050239	48	46

Plot of number of reads by sample

sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)

p <- ggplot(metadata_read, aes(x = Sample, y = Read_depth))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90)) +
  ggtitle("") + 
  theme(legend.title = element_text(size = 18), 
        legend.position="bottom",
        legend.text=element_text(size=11), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        strip.text.x = element_text(size = 12)) +
 facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=4, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 550000)+
  geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=2, hjust=-0.25, vjust=0.25, angle=90)+
  guides(fill=guide_legend(title="Oligotype", label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))

p

Filter samples

Remove sample with number of reads <1000

a <- prune_samples(sample_sums(ps.filter)<=1000, ps.filter)
test <- data.frame(a@sam_data)
test <- test[test$Sample!="NP20" & test$Sample!="NP29" & test$Sample!="NP30" & test$Sample!="NP34" & test$Sample!="NP36",]
toremove <- test$Sample

ps.filter <- subset_samples(ps.filter, !(Sample %in% toremove))
ps.filter

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 67 taxa and 113 samples ]
## sample_data() Sample Data:       [ 113 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 67 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 67 reference sequences ]

compute_read_counts(ps.filter)

## [1] 6452623

Plot of final number of reads by sample

sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)
metadata_read_whole <- metadata_read[metadata_read$Organ=="Whole",]
metadata_read_ovary <- metadata_read[metadata_read$Organ=="Ovary",]


guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))

p <- ggplot(metadata_read_whole, aes(x = Sample, y = Read_depth))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 18), 
        legend.position="bottom",
        legend.text=element_text(size=16), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        strip.text.x = element_text(size = 16),
        axis.text = element_text(size = 15)) +
 facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 1000000)+
  geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
  guides(fill=guide_legend(title="Oligotype"))

p

p2 <- ggplot(metadata_read_ovary, aes(x = Sample, y = Read_depth))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 18), 
        legend.position="bottom",
        legend.text=element_text(size=16), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        strip.text.x = element_text(size = 16),
        axis.text = element_text(size = 15)) +
 facet_wrap(~Species_italic+Strain_italic+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 1000000)+
  geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
  guides(fill=guide_legend(title="Oligotype"))

p2

# panels
p_group <- plot_grid(p+theme(legend.position="none"), 
          p2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A", "B"), c(0, 0), c(1, 0.5), size = 15)

legend_plot <- get_legend(p + theme(legend.position="bottom"))

p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Counts after removing samples

Create a dataframe

x <- c("Culex pipiens", "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar", 
       "Culex quinquefasciatus", "Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)",
       "Aedes aegypti (Guadeloupe)",
       "Total")
y <- c("Reads", "Oligotypes", "Samples")

df <- data.frame(matrix(ncol=3, nrow=9))
rownames(df) <- x
colnames(df) <- y

df2 <- df
df3 <- df

Whole + Ovary

# Culex pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)

## All Strain
### oligotype
nrow(ps.pipiens@otu_table) -> df["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data) -> df["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens) -> df["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Bosc")) -> df["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Field - Camping Europe")) -> df["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Strain=="Laboratory - Lavar")) -> df["Laboratory - Lavar", "Reads"]


# Culex quinquefasciatus
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)

## All Strain
### oligotype
nrow(ps.quinque@otu_table) -> df["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data) -> df["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque) -> df["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque %>% 
                 subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
                 check_ps())) -> df["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df["Laboratory - Slab TC (Wolbachia -)", "Reads"]


# Aedes aegypti
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes %>% subset_samples(Strain=="Field - Guadeloupe")) -> df["Aedes aegypti (Guadeloupe)", "Reads"]

# Total
### oligotype
nrow(ps.filter@otu_table) -> df["Total", "Oligotypes"]
### samples
nrow(ps.filter@sam_data) -> df["Total", "Samples"]
### reads
compute_read_counts(ps.filter) -> df["Total", "Reads"]

df %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	2574050	49	82
Field - Bosc	803574	36	23
Field - Camping Europe	869064	38	24
Laboratory - Lavar	901412	48	35
Culex quinquefasciatus	1924657	59	20
Field - Guadeloupe	1721053	47	7
Laboratory - Slab TC (Wolbachia -)	203604	45	13
Aedes aegypti (Guadeloupe)	1953916	54	11
Total	6452623	67	113

Whole

# Culex pipiens
ps.pipiens.whole <- ps.pipiens %>% subset_samples(Organ=="Whole")
ps.pipiens.whole <- ps.pipiens.whole %>% check_ps()

## All Strain
### oligotype
nrow(ps.pipiens.whole@otu_table) -> df2["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data) -> df2["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole) -> df2["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df2["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Bosc")) -> df2["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df2["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Field - Camping Europe")) -> df2["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.whole %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df2["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Strain=="Laboratory - Lavar")) -> df2["Laboratory - Lavar", "Reads"]



# Culex quinquefasciatus
ps.quinque.whole <- subset_samples(ps.quinque, Organ=="Whole")
ps.quinque.whole <- check_ps(ps.quinque.whole)

## All Strain
### oligotype
nrow(ps.quinque.whole@otu_table) -> df2["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data) -> df2["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.whole) -> df2["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.whole %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df2["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque.whole %>% 
                 subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
                 check_ps())) -> df2["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df2["Laboratory - Slab TC (Wolbachia -)", "Reads"]



# Aedes aegypti
ps.aedes.whole <- subset_samples(ps.aedes, Organ=="Whole")
ps.aedes.whole <- check_ps(ps.aedes.whole)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.whole %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df2["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.whole@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.whole %>% subset_samples(Strain=="Field - Guadeloupe")) -> df2["Aedes aegypti (Guadeloupe)", "Reads"]


# Total
ps.whole <- ps.filter %>% subset_samples(Organ=="Whole")
ps.whole <- ps.whole %>% check_ps()
### oligotype
nrow(ps.whole@otu_table) -> df2["Total", "Oligotypes"]
### samples
nrow(ps.whole@sam_data) -> df2["Total", "Samples"]
### reads
compute_read_counts(ps.whole) -> df2["Total", "Reads"]

df2 %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	1251213	45	41
Field - Bosc	545790	33	13
Field - Camping Europe	98251	33	6
Laboratory - Lavar	607172	45	22
Culex quinquefasciatus	205364	57	17
Field - Guadeloupe	1760	38	4
Laboratory - Slab TC (Wolbachia -)	203604	45	13
Aedes aegypti (Guadeloupe)	1945807	54	9
Total	3402384	67	67

Ovary

# Culex pipiens
ps.pipiens.ovary <- ps.pipiens %>% subset_samples(Organ=="Ovary")
ps.pipiens.ovary <- ps.pipiens.ovary %>% check_ps()

## All Strain
### oligotype
nrow(ps.pipiens.ovary@otu_table) -> df3["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data) -> df3["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary) -> df3["Culex pipiens", "Reads"]


## Bosc
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Field - Bosc") %>% 
                 check_ps())) -> df3["Field - Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Bosc")) -> df3["Field - Bosc", "Reads"]


## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Field - Camping Europe") %>% 
                 check_ps())) -> df3["Field - Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Field - Camping Europe")) -> df3["Field - Camping Europe", "Reads"]

## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.ovary %>% 
                 subset_samples(Strain=="Laboratory - Lavar") %>% 
                 check_ps())) -> df3["Laboratory - Lavar", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Strain=="Laboratory - Lavar")) -> df3["Laboratory - Lavar", "Reads"]



# Culex quinquefasciatus
ps.quinque.ovary <- subset_samples(ps.quinque, Organ=="Ovary")
ps.quinque.ovary <- check_ps(ps.quinque.ovary)

## All Strain
### oligotype
nrow(ps.quinque.ovary@otu_table) -> df3["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data) -> df3["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary) -> df3["Culex quinquefasciatus", "Reads"]


## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.ovary %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df3["Field - Guadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Field - Guadeloupe", "Reads"]


## Wolbachia -
### oligotype
# nrow(otu_table(ps.quinque.ovary %>% 
#                  subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)") %>% 
#                  check_ps())) -> df3["Laboratory - Slab TC (Wolbachia -)", "Oligotypes"]
### samples
# nrow(ps.quinque.ovary@sam_data %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Samples"]
# ### reads
# compute_read_counts(ps.quinque.ovary %>% subset_samples(Strain=="Laboratory - Slab TC (Wolbachia -)")) -> df3["Laboratory - Slab TC (Wolbachia -)", "Reads"]



# Aedes aegypti
ps.aedes.ovary <- subset_samples(ps.aedes, Organ=="Ovary")
ps.aedes.ovary <- check_ps(ps.aedes.ovary)

## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.ovary %>% 
                 subset_samples(Strain=="Field - Guadeloupe") %>% 
                 check_ps())) -> df3["Aedes aegypti (Guadeloupe)", "Oligotypes"]
### samples
nrow(ps.aedes.ovary@sam_data %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Samples"]
### reads
compute_read_counts(ps.aedes.ovary %>% subset_samples(Strain=="Field - Guadeloupe")) -> df3["Aedes aegypti (Guadeloupe)", "Reads"]


# Total
ps.ovary <- ps.filter %>% subset_samples(Organ=="Ovary")
ps.ovary <- ps.ovary %>% check_ps()
### oligotype
nrow(ps.ovary@otu_table) -> df3["Total", "Oligotypes"]
### samples
nrow(ps.ovary@sam_data) -> df3["Total", "Samples"]
### reads
compute_read_counts(ps.ovary) -> df3["Total", "Reads"]

df3[is.na(df3)] <- 0

df3 %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

	Reads	Oligotypes	Samples
Culex pipiens	1322837	34	41
Field - Bosc	257784	26	10
Field - Camping Europe	770813	27	18
Laboratory - Lavar	294240	32	13
Culex quinquefasciatus	1719293	26	3
Field - Guadeloupe	1719293	26	3
Laboratory - Slab TC (Wolbachia -)	0	0	0
Aedes aegypti (Guadeloupe)	8109	30	2
Total	3050239	48	46

Save data

Save count tables and count plots

write.table(df, "../../../output/1_MED/1D/1D_Counts_all.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(df2, "../../../output/1_MED/1D/1D_Counts_whole.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(df3, "../../../output/1_MED/1D/1D_Counts_ovary.tsv", sep="\t", row.names = TRUE, col.names=NA)

tiff("../../../output/1_MED/1D/1D_counts_by_sample.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()

## quartz_off_screen 
##                 2

png("../../../output/1_MED/1D/1D_counts_by_sample.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()

## quartz_off_screen 
##                 2

Save metadata (for Supplementary Table 1)

count_after_decontam <- data.frame(ps.filter@otu_table)
tax_after_decontam <- data.frame(ps.filter@tax_table)
metadata_after_decontam <- data.frame(ps.filter@sam_data)

compute_read_counts(ps.filter)

## [1] 6452623

write.table(count_after_decontam, "../../../output/1_MED/1D/1D_count_supplementary_table1.tsv", sep="\t", row.names = TRUE, col.names=NA)
write.table(tax_after_decontam, "../../../output/1_MED/1D/1D_tax_supplementary_table1.tsv", sep="\t", row.names = TRUE, col.names=NA)

Save phyloseq object after decontam

saveRDS(ps.filter, "../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")

Session info

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.0    microbiome_1.8.0 kableExtra_1.3.1 decontam_1.6.0  
##  [5] forcats_0.5.0    stringr_1.4.0    dplyr_1.0.2      purrr_0.3.4     
##  [9] readr_1.4.0      tidyr_1.1.2      tibble_3.0.4     tidyverse_1.3.0 
## [13] ggplot2_3.3.2    phyloseq_1.30.0 
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-149        fs_1.5.0            lubridate_1.7.9    
##  [4] webshot_0.5.2       httr_1.4.2          tools_3.6.3        
##  [7] backports_1.1.10    bslib_0.2.4         R6_2.4.1           
## [10] vegan_2.5-6         DBI_1.1.0           BiocGenerics_0.32.0
## [13] mgcv_1.8-33         colorspace_2.0-0    permute_0.9-5      
## [16] ade4_1.7-15         withr_2.3.0         tidyselect_1.1.0   
## [19] compiler_3.6.3      cli_2.1.0           rvest_0.3.6        
## [22] Biobase_2.46.0      xml2_1.3.2          labeling_0.4.2     
## [25] bookdown_0.22       sass_0.3.1          scales_1.1.1       
## [28] digest_0.6.26       rmarkdown_2.7       XVector_0.26.0     
## [31] pkgconfig_2.0.3     htmltools_0.5.1.1   highr_0.8          
## [34] dbplyr_1.4.4        rlang_0.4.10        readxl_1.3.1       
## [37] rstudioapi_0.11     farver_2.0.3        jquerylib_0.1.3    
## [40] generics_0.0.2      jsonlite_1.7.1      magrittr_1.5       
## [43] biomformat_1.14.0   Matrix_1.2-18       fansi_0.4.1        
## [46] Rcpp_1.0.5          munsell_0.5.0       S4Vectors_0.24.4   
## [49] Rhdf5lib_1.8.0      ape_5.4-1           lifecycle_0.2.0    
## [52] stringi_1.5.3       yaml_2.2.1          MASS_7.3-53        
## [55] zlibbioc_1.32.0     Rtsne_0.15          rhdf5_2.30.1       
## [58] plyr_1.8.6          grid_3.6.3          blob_1.2.1         
## [61] parallel_3.6.3      crayon_1.3.4        lattice_0.20-41    
## [64] Biostrings_2.54.0   haven_2.3.1         splines_3.6.3      
## [67] multtest_2.42.0     hms_0.5.3           ps_1.4.0           
## [70] knitr_1.30          pillar_1.4.6        igraph_1.2.6       
## [73] reshape2_1.4.4      codetools_0.2-16    stats4_3.6.3       
## [76] reprex_0.3.0        glue_1.4.2          evaluate_0.14      
## [79] data.table_1.13.2   modelr_0.1.8        vctrs_0.3.4        
## [82] rmdformats_1.0.2    foreach_1.5.1       cellranger_1.1.0   
## [85] gtable_0.3.0        assertthat_0.2.1    xfun_0.22          
## [88] broom_0.7.2         viridisLite_0.3.0   survival_3.2-7     
## [91] iterators_1.0.13    IRanges_2.20.2      cluster_2.1.0      
## [94] ellipsis_0.3.1